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INTRODUCTION 

This  report  is  a continuation  of  Report  Number  13, 
"Parametric  Optimal  Design  Part  I:  First  Order  Algorithm". 

In  this  report,  two  new  methods  of  parametric  optimal  design 
(POD)  are  developed.  The  first  method  is  called  a second 
order  algorithm,  and  the  second  one  is  called  a hybrid 
algorithm.  The  hybrid  algorithm  is  a combination  of  the 
first  and  the  second  order  methods.  Chapter  ^ presents  a 
design  sensitivity  analysis  with  error  compensation  which 
makes  second  order  algorithms  possible.  Hybrid  methods  are 
discussed  in  Chapter  5-  Several  example  problems  are  solved, 
using  both  the  second  order  algorithm  and  the  hybrid  method, 
and  a comparison  of  results  is  made.  In  Chapter  6,  detailed 
procedures  for  implementation  of  various  algorithms  is  pre- 
sented. Chapter  7 provides  a summary  of  the  main  contribu- 
tions and  comparison  of  various  methods  alongwith  the 
difficulties  involved  in  applications  of  the  developed 
algorithms . 
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CHAPTER  4 

SECOND  ORDER  ALGORITHM 
4.1  Introduction 

The  first  order  algorithm  described  in  the  previous 
chapter  is  equivalent  to  alternate  maximization  and  mini- 
mization. This  is  a feasible  method  and  has  been  used 
in  solving  the  classical  minimax  problem.  It  is  apparent, 
however,  that  during  the  maximization  procedure  in  the 
subproblem,  if  one  can  account  for  higher  order  changes 
in  the  design  variables,  one  will  have  a better  estimate 
of  a for  the  next  iteration,  without  a separate  maximi- 
zation procedure.  In  other  words,  one  wishes  to  use  a 
higher  order  approximation  to  predict  the  effect  on  the 
maximum  point  of  the  subproblem,  considering  the  change 
in  the  design  variable  6b.  Furthermore,  if  the  current 
maximum  point  a is  not  exact,  then  this  error  should  be 
corrected  In  the  next  iteration.  The  algorithm  developed 
in  this  section  is  based  on  this  Idea. 

The  procedure  to  be  employed  is  divided  into  two  major 
steps,  (1)  Sensitivity  analysis  for  the  subproblem  and 
(2)  a minimization  procedure  for  the  outer  problem,  using 
sensitivity  information. 
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^.2  Sensitivity  Analysis  with  Error  Compensation 

Conventionally,  sensitivity  analysis  consists  of  de- 
termination of  sensitivity  of  the  objective  function  with 
respect  to  a change  in  constraint  level,  and  determination 
of  the  sensitivity  of  the  solution  point  to  a change  in 
the  constraint. 

A more  general  sensitivity  analysis  includes  the 
sensitivity  of  the  solution  point  with  respect  to  a set  of 
parameters  that  appear  in  the  problem.  This  last  type  is 
of  concern  for  POD,  and  will  be  referred  to  as  sensitivity 
analysis  unless  otherwise  specified.  As  explained  in  the 
introduction  to  this  chapter,  if  one  obtains  a functional 
relation  for  the  solution  point  of  subproblem,  a * a(b), 
the  POD  problem  can  be  reduced  to  a classical  nonlinear 
programming  problem  by  substitution  of  this  relation  for 
each  subproblem.  This  relation  may  be  considered  a closed 
loop  solution  of  the  subproblem,  or  a feedback  law  for  the 
maximum  point.  Theoretically,  then,  the  POD  problem  can 
be  solved  by  obtaining  the  closed  loop  solution.  However, 
the  closed  loop  solution  can  not  be  found  in  most  practical 
problems . 

The  approach  that  is  persued  here  involves  solution 
for  a linear  approximation  to  the  feedback  law,  in  the 
neighborhood  of  the  current  solution  point.  The  sensitiv- 
ity of  the  maximum  point  a,  with  respect  to  a change  in  b. 
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is  then  Riven  by  the  first  order  feedback  law  and  is  ex- 
pressed by  — 

Sensitivity  analysis  for  a general  nonlinear  program- 
ming problem  is  considered  by  Fiacco  [42],  The  problem 
treated  is  to  find  the  partial  derivatives  of  the  coordi- 
nates of  the  solution  point  and  the  objective  function 
with  respect  to  the  free  parameters.  These  derivatives 
are  called  the  first  order  sensitivity  of  a second  order 
local  solution.  It  is  shown  in  [42]  that  under  the  as- 
sumptions of  a certlan  degree  of  differentiability  and 
some  regularity  conditions,  the  derivatives  are  well  de- 
fined . 

The  results  of  [42]  might  be  used  in  obtaining  explic- 
it sensitivity  information  for  the  problem  at  hand,  pro- 
vided the  solution  point  and  the  Lagrange  multipliers  are 
known.  However,  when  there  are  many  state  equations,  this 
sensitivity  analysis  requires  excessive  computer  storage 
and  computer  time.  Therefore,  in  subsequent  papers  [43,44] 
the  same  authors  use  an  algorithm  of  direct  calculation  of 
the  components  of  the  sensitivity  derivatives,  using  finite 
differences . 

In  this  section,  an  explicit  expression  for  the  first 
order  feedback  law  is  obtained  by  using  a gradient  projec- 
tion method  with  constraint  error  compensation.  This  meth- 
od is  derived  by  use  of  a second  order  approximation  of  the 


objective  function  of  the  subproblem  Eq.  (2.3)  and  of  the 
state  equations  (2.2).  A first  order  approximation  of  the 
a-constraint  equations  (2.5)  is  employed.  The  second  or- 
der approximation  is  essential,  since  with  a linear  ap- 
proximation, the  maximum  point  a is  independent  of  b and 
the  sensitivity  is  zero. 

The  first-  order  feedback  law  is  now  obtained  by  devel 
oping  an  explicit  feedback  law  for  the  second  order  ap- 
proximate problem. 

The  inner  problem  is  of  the  following  type:  Choose  a 

to  maximize 


g(z ,b  ,a ) 

(4.1) 

subject 

to 

h(z,b,a)  = 0, 

(4.2) 

q (a  ) <_  0 . 

(4.3) 

Since  z 

is  considered  as  a function  of  b and  a , 

one  can 

express  the  variation  of  g in  terms  of  variations  of  b and 
a , as  follows ; 

T T IT  T 

fig  * B fib  + A 5a  + ^fib  BBfib  + fib  BAfia 

+ ^fia*^  AAfia  . (4.4) 


Also , 

5 q = *—  6 a = Q^fi  a < A Q , (4.5) 

3a 

where  ( ) denoces  the  constraints  that  are  tight  or  violat- 
ed and  Aq  is  the  desired  reduction  in  constraint  violation 
usually  taken  Aq=  -q(a").  The  various  matrices  in  the 
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expression  for  6g  are  given  in  Appendix  D.  In  Eq.  (4.5), 
only  a linear  approximation  of  the  a constraints  has  been 
used.  In  many  problems,  q(a)  is  linear  in  a. 

The  approximate  problem  is  now  to  maximize  fig  of  Eq. 

(4.4)  with  respect  to  6a,  subject  to  the  constraint  of  Eq. 

(4.5) .  To  Insure  the  validity  of  the  approximations,  a 
restriction  is  placed  on  6a  such  that 


where  WQ  is  a positive  definite  weighting  matrix,  and  is 

predetermined  in  the  same  way  as  explained  for  W in  the 

b 

previous  chapter. 

Following  Section  2.4.4,  the  Kuhn-Tucker  necessary 
conditions  for  the  problem  (4.4),  (4.5),  and  (4.6),  where 
6a  is  now  the  design  variable  as  an  NLP,  are  solved  to 
obtain , 


(4.6) 


a 


6a  = -D6b  - C 


(4.7) 


where 


D = ♦ _1{ I-Mqq } • BAT , 

C = ♦ -1[ { I-M  ) • A - M Aq] 


qq  q 


(4.8) 


(4.9) 


and 


(4.10) 


where 


o = AA  - 2vW 

i 

M = QM'1 , 

Q -t  -1- 

M = Q1#  Q, 


a 


(4.11) 


(4.12) 
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~T  -1 
M = M Q ♦ . 

qq  q 

The  constants  v and  p are  multipliers  to  Eqs.  (4.6)  and 
(4.5),  respectively. 

Eq . (4.7)  is  called  the  first  order  feedback  law  with 

error  compensation,  where  -C  is  the  compensating  term  for 

the  error  in  the  nominal  solution  to  the  inner  problem. 

A necessary  condition  for  the  nominal  solution  to  be  the 

solution  of  the  subproblem  is  that  C=0  and  $ is  negative 

semi-definite,  since  ♦ is  the  Hessian  of  the  Lagrangian 

functional  to  problem  (4.4)  and  (4.5)  (Section  2.4.2). 

Hence,  in  a neighborhood  of  the  maximum  point,  the  feedback 

8 (x 

law  is  given  by  6a»  -D6b  and  the  sensitivity  is  given 
by  -D. 

The  multiplier  v is  determined  as  a stepsize,  instead 
of  calculating  it  as  the  Lagrange  multiplier  for  the  step- 
size  restriction  of  Eq.  (4.6).  This  approach  is  employed, 
since  the  stepsize  restriction  contains  another  constant 
t that  has  little  more  physical  significance  than  v itself. 
The  choice  of  v will  be  discussed  in  a later  section.  By 
choosing  v>0,  ♦ can  be  made  negative  definite  and  inversion 
of  ♦ is  always  feasible.  This  will  be  recognized  as  a 
great  advantage  to  numerical  computation,  even  for  the  case 
in  which  AA  is  nearly  singular.  For  a system  without  con- 
straint violation,  0*”^A  is  simply  the  correction  obtained 
by  Newton’s  method,  with  a stepsize  restriction.  The 


87 


weighting  matrix  WQ  changes  the  local  scales  in  design 
space.  A wise  choice  of  this  matrix  can  improve  the  con- 
vergence of  the  overall  iterative  method.  These  points 
are  discussed  in  Chapter  6 when  the  Newton's  method  is  de- 


a-space  onto  the  tangent  subspace  determined  by  the  q 
constraint  set.  Thus,  C is  a vector  from  A,  projected 
onto  the  tangent  subspace  and  then  adjusted  by  a local 
metric  Wc . Similarly,  the  sensitivity  D is  the  matrix 


transformed  from  the  cross  derivative  BA  by  the  projection 
matrix  and  the  local  metric  W„ . 

In  the  follwoing  section,  Eq.  (4.7)  is  utilized  to 
obtain  a second  order  approximation  of  the  POD  problem,  in 
terms  of  fib,  and  to  solve  the  outer  problem. 

4.3  Solution  of  the  Second  Order  Approximate  Problem 
The  minimization  of  6f,  subject  to  the  reduced  con- 
straints, without  the  environmental  parameter,  will  be 
termed  the  outer  procedure  for  the  POD  problem.  By  direct 
substitution  from  the  preceding  section,  the  POD  problem 
can  be  approximated  explicitly  in  terms  of  fib  as:  Minimize 


rived . 


It  can  be  shown  by  direct  multiplication  [34]  that 
(I-M  ) is  a projection  matrix  and  that  it  projects  the 


-T 


6 f ■ fib  + |-5bTLJfib, 


(4.13) 


subject  to 


(4.14) 
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where 


a T m . J 3 f 

tj  » (B-DiA-BA*C+Di -AA-Oj  , l = 

* [gB+DT- AA-D-(BA-D+DT*BAT)] . , LJ  » 
J J 


(4.15) 


K j = (ATC-  |cT- AA-Oj,  J e Ie  , 


and  Ag=  -g  is  the  intended  reduction  in  the  current  con- 
straint violation.  The  superscript  <f>  denotes  the  inclu- 
sion of  all  the  e -active  constraints  whose  indices  are  in 
I£  . It  is  noted  that  the  vector  in  Eq . (4.14)  is  a 
formal  expression  to  the  total  derivative  of  g with  respect 
to  b,  when  C is  assumed  zero.  If  the  optimality  condition 
for  the  subproblem  is  utilized,  the  expression  for  is 
simplified  and  K*  is  zero,  since  C=  0 in  this  case.  But  for 
numerical  calculation,  this  optimality  condition  is  satis- 
fied only  approximately.  Hence,  it  is  necessary  to  retain 
C. 

To  solve  the  second  order  approximate  problem,  the 
Newton-Raphson  technique  is  employed  to  solve  the  Kuhn- 
Tucker  necessary  conditions.  To  assure  that  the  approxi- 
mate problem  is  close  to  the  original  problem,  an  additional 
stepsize  restriction  is  imposed  on  6b.  Prom  Eq.  (4.7) 


I | I*"1!  M I I-Mqq|  | • | | BAT|  | -|  |6b|  | + | | C|  where 


| | $ | | denotes  a matrix  norm,  defined  by  ||»||a  ||$x||'. 


and  | | • | | and  | | • | | ' are  norms  on  Rn  and  Rp . For 
a normal  problem,  $ is  negative  definite  and  the  spectral 


=1 
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norm  of  the  projection  matrix 
to  1.  Therefore, 


(I-M  ) 

qq 


is  less  than  or  equal 


iu«ir  i Kii«bii  + 1 ici  r,  (n.18) 

where  K= | | o-1 | | • | | I-M  ||*||bAT||  is  finite.  Since 
llcir-0  as  the  nominal  solution  approaches  the  solution 
to  the  subproblem,  if  | 1 6 b | | is  small,  then  '|6a||f  is  small; 
but  not  necessarily  vice  versa.  Hence,  the  stepsize  re- 
striction on  | 1 6a  | |',  given  by  Eq.  (4.6),  is  not  sufficient 
and  the  following  stepsize  restriction  is  imposed  on  6b: 

6bT  W,  6b  < n2  ( ^ . IS ) 

D “ 

where  is  a positive  definite  design  variable  weighting 
matrix,  predetermined  as  described  in  Section  6.3.2. 

The  Kuhn-Tucker  conditions  for  the  approximate  problem 
are,  using  summation  notation, 

li  + Lij{bJ  + 'tj\f  + L}jk«bj»k  * SrWbijtbj  = o, 

(4.20) 

X1(£j16bj+  ^b^^ebj  - k[  - Ag^  = 0,  (4.21) 

y(6b^W^i^6bj  - n2)  * 0,  (4.22) 

where  X^  >_  0,  and  y >_  0 are  multipliers. 

The  system  of  Eqs.  (4.20),  (4.21),  and  (4.22)  is  high- 
ly nonlinear  in  the  unknowns  6b,  X,  and  y.  The  existence 
of  a solution  is  difficult  to  prove.  Even  when  a solution 
exists  for  a physically  well  formulated  problem,  it  may  not 
be  unique.  Further,  because  of  the  sign  restriction  on 
multipliers,  the  proper  solution  is  difficult  to  obtain. 


in  practice.  If  a solution  of  the  system  of  equations 
yields  a negative  multiplier,  the  corresponding  constraint 
is  discarded  and  the  computation  is  repeated.  To  solve  the 
resulting  system,  the  Newton-Raphson  method  is  used,  with 
the  following  derivative  matrix: 


LiJ+xkLiJk+2YWb1J » 1ij+LikJ6bk 


Xi(£ji+Lkji6bk) 


’ 6ijFn+l 


(i  not  summed) 

2Y6biWblj  , 0 


2Wh  6b, 
biJ  J 


6b^Wb6b-n^  J 


(4.23) 


where 

Fn+i  * *jiabJ  + \ 6bkLkJi6bJ  - Ki  - A*i* 

This  derivative  matrix  is  nonsymmetric*  and  nonsingular, 

as  long  as  the  system  is  normal. 

The  success  of  the  Newton-Raphson  method  depends 
largely  on  the  initial  estimate.  Since  the  solution  for 
6b  is  expected  to  be  near  the  origin  of  design  space, 

6b  3 0 is  found  to  be  effective  as  the  initial  estimate. 

The  initial  choice  of  X and  y i3  obtained  by  the  algo- 
rithm given  in  Section  2.4.4,  with  and  interpreted 

A 

as  given  in  the  present  section  and  Ag  replaced  by  Ag+KT . 


+If  equality  is  presumed  to  hold  in  Eqs.  (4. 14),  and 
(4.19)  instead  of  Eqs.  (4.21),  and  (4.22),  the  derivative 
matrix  would  be  symmetric.  But  the  solution  set  from  this 
system  of  equations  may  be  too  restrictive.  It  is  only  a 
subset  of  the  solution  set  of  the  original  system,  (4.20), 
(4.21),  and  (4.22). 
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The  vectors  6b  and  x are  given  as: 

.1 


6b 


6b' 


2Y 

-1 


+ 6 b 


W 


= - 27~(t  nV  )+  w;Vx2, 


(1* . 24) 


where 


X1  - 2Yl2 


(t*  W^tMx1  = -l* 
U*  W^t^X2 


(4.25) 

(4.26) 


: Kr  + Ag. 

It  is  noted  that  this  Is  different  from  the  result  of  the 
pure  first  order  algorithm,  as  given  in  Section  3.2,  since 

the  expression  of  £*  in  this  section  is  different  from 

4 T 

£ -{B.,  Je I } in  the  first  order  algorithm.  Also,  in 
Section  3.2,  K*  =0. 

It  is  shown  in  Appendix  E that  the  Lagrange  multiplier 
for  the  first  order  approximate  problem  is  the  same  as  for 
the  original  problem.  The  Lagrange  multiplier  for  the  sec- 


ond order  approximate  problem  is  obtained  as 

T T 

ql1  = +5b'TL*  J+LJ  fb), 

Q>2  * -itf  <K*+*g>, 


where 


M 


♦ ♦ 


T T 

(t*  +6bTL+  )VT1(£*+L*  6¥), 
b 


F 


6 kI 


5TT 


(4.27) 

(4.28) 

(4.29) 

(4.30) 


and  , if  , £^,  and  L^  are  claculated  at  the  solution  of 
the  approximate  problem,  6"b.  Comparing  with  the  above 
first  order  result,  Eq . (4.25),  qX^  is  equal  to  X1  when 
6b  = 0.  In  the  limit,  the  solution  of  the  second  order 
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approximate  problem  Is  zero.  During  iteration,  however, 
fib  is  nonzero  and  ^A^/A^.  Therefore,  if  the  ini- 

tial estimate  for  the  Newton-Raphson  method,  5b^°^,  is 
zero,  A and  A from  Eqs.  (4.25)  and  (4.26)  are  the  proper 
estimates  for  ^x.  As  each  outer  porcedure  converges,  the 
solution  of  the  second  order  approximate  problem  must  con- 
verge to  zero.  As  the  iterate  of  POD  converges,  ^x^ax1 

q 2 2 

and  A = A , since  fib -*0 . Further,  as  shown  in  Appendix  h, 

1 2 

A and  A‘" , from  Eqs  (4.25)  and  (4.26),  approach  the  Lagrange 
multipliers  of  the  POD  problem. 

4.4  Numerical  Examples 

The  weapons  allocation  problem  described  in  Section 
2.3.1  is  solved  by  the  second  order  method  of  this  chapter. 
With  the  same  initial  estimates  as  given  in  Section  3.3, 
the  solution  is  obtained  only  after  8 iterations  (Table  7); 
b * [0.0,  0.18186,  0.27273]T,  with  a cost  of  -2.2340  and 

m 

the  worst  case  a = [0.0,  0.1751,  0.4654]  . As  compared 
with  the  behaviour  of  iterates  in  the  first  order  method 
(Table  1),  Table  7 shows  that  convergence  is  rapid  as  the 
Iterates  approach  the  solution.  Since  the  problem  is  con- 
vex, a quadratic  rate  of  convergence  can  be  expected.  The 
other  feature  is  that  the  maximization  procedure  for  the  in- 
ner problem,  which  is  an  essential  part  of  the  first  order 
algorithm,  is  not  necessary;  since  the  second  order  method 
accounts  for  the  change  in  maximum  points,  as  well  as 
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compensating  any  current  error  in  the  nominal  solution  of 
inner  problems. 

The  sensitivity  coefficient,  -D  in  Eq . (4.7),  at  the 
solution,  is  found  to  be 

' 3. 7 3 ( — 8 ) 2.61  (-8)  1.86(-8) 

9 . 06(  -2)  4.55(-l)  7 . 70(  -2) 

9.7M-?)  7 - 7 0 ( — 2 ) 3 . 5 8 ( —1 ) _ 

Tables  8 and  9 show  the  solutions  of  the  vibration 
isolator  design  problem  (Section  2.3.3),  with  frequency 
ranges  of  [0.8,  1.5]  and  [0.9,  1.25],  respectively.  The 
initial  estimate  for  case  (i)  is  rather  good  and  the  method 
converges  rapidly.  On  the  other  hand,  the  initial  estimate 
for  case  (ii)  is  rather  bad.  Since  the  artificial  design 
variable  has  been  used,  the  first  6 iterations  only  adjust 
the  artificial  design  variable.  After  Iteration  number  6, 
rather  good  convergence  is  obtained.  It  is  noted  that 
case  (ii)  behaves  badly,  with  regard  to  the  system  of  Kuhn- 
Tucker  equations,  (4.20),  (4.21),  and  (4.22).  After  fix- 
ing y>0,  determined  by  the  method  described  in  Section 
6.3*2,  the  Newton’s  method  converges  very  rapidly.  This 
may  limit  the  convergence  rate  of  outer  procedure,  but  the 
precise  affect  is  not  known.  The  computer  time  per  itera- 
tion on  an  IBM  360/65  was  about  0.24  seconds,  in  FORTRAN(H). 

Sensitivity  coefficients  are  obtained  as:  For  case 

(1),  -D  = [0.835,  0.197]  at  a=0.848,  and  -D  = 
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[-1.103,  0.300]  at  a=1.059.  For  case  (11)  -D  = 

C 1 - 36 ( —9 ) , -8. 05(-8)]  at  o = 0.9,  and  -D  = [-0.574,  0.657] 
at  a=l.H9. 

McMunn's  problem  [32],  described  In  Section  2.3*5,  Is 
solved,  with  results  In  Table  10.  This  problem  has  a 
saddle  point  solution.  As  pointed  out  earlier,  however, 
since  one  is  not  sure,  at  the  beginning,  how  many  maxima 
there  are  for  a given  range  of  free  parameters,  a certain 
amount  of  Intuition  is  necessary.  Four  initial  estimates 
of  maximum  points  are  used  here.  Table  10  shows  that  they 
merge  into  one  point,  meaning  one  maximum  point  at  the 
optimum  design.  The  same  initial  design  estimates,  as 
given  in  [32],  are  used.  As  in  the  vibration  isolator  de- 
sign case  (ii)  (Table  9),  the  first  7 iterations  are  spent 
compensating  the  constraint  violations.  After  the  7-th 
iteration,  it  is  seen  that  the  rate  of  convergence  is 
quadratic.  The  final  solution  is  b = [ 2 . 099  ,0 . 9133 ,1 • 145] * , 
with  J = 1.49-0  attained  at  a=1.505.  The  computer  time  was 
about  1 second  per  iteration.  McMunn's  solution  was 
b = [2. 1160,0. 9149, 1.0695]T,  and  J = 1.494  at  a=1.503. 

Note  that  part  of  the  difference  may  come  from  the  round- 
off error  in  the  numerical  values  of  the  matrices  M,  D,  and 
K.  A sensitivity  coefficient  of  -D  = [0.976,5-414,0.4^7] 
is  obtained  during  the  solution  procedure. 
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In  the  notation  of  Section  2.2,2,  the  following 
Chebyshev  approximation  problem  is  taken  from  [19]; 


ru>=  ii.  • 

and 

b^a  bha 

F(b,a)  = b^e  + b£e  , a c [0.,  1.0]. 

Table  11 

shows  the  results.  In  [19],  the  solution  is 

obtained 

as  b={0. 713835,  0.285966,  -0.407327,  -2.442953), 

with  the 

square  of  maximum  error  in  between  0.23(-7)  and 

0 . 76 ( -7 ) . 

There  is  some  difference  in  the  solution, 

although 

||<5b^||  Is  quite  small,  satisfying  any  practical 

stopping 

criterion.  It  was  found  that  even  though  the 

iterations  were  continued  no  significant  improvement  was 
obtained.  This  may  be  caused  by  round-off  errors.  The 
computer  time  was  about  0.33  seconds  per  iteration. 


Table  7 Weapons  Allocation  Problem 
(Second  Order  Algorithm) 


iter 

bi 

obj  . 

1 1 Sb1 | | 

ai 

, tight 

i 1 6a  ||  Si,  q. 

0.1 

0.1 

0.2 

0.3 

0 

0.1 

0.2 

0.3 

-2.10 

0.795 

0.0 

0.246 

0 . 563 

(4)0.002  q 

1 

-0.579 

0.323 

0.481 

-2.48 

0.23 

0.0 

0.219 

0.533 

62,ql 

2 

0.0 

0.309 

0.468 

-2.014 

0.91 

0.0 

0.219 

0.533 

(0)0.15 

3 

0.0 

0.193 

0.277 

-2.233 

1.0 

0.0 

0.193 

0.486 

(2)0.004 

4 

0.0 

0.193 

0.277 

-2.233 

0.14 

0.0 

0.194 

0.482 

(0)0.003 

c; 

0.0 

0.178 

0.270 

-2.235 

0.02 

0.0 

0.183 

0.473 

(0)0.008 

6 

0.0 

0.181 

0.271 

-2.234 

0.01 

0.0 

0.176 

0.467 

(0)0.001 

7 

0.0 

0.182 

0.272 

-2.234 

0.003 

0.0 

0.175 

0.465 

(1)0.0002 

8 

0.0 

0.18186 

0.27273 

-2.2340 

0.0002 

0.0 

0.1751 

0.4654 

1 

V-/ 

ro 

0 

O 

V_/ 

VO 

o\ 


Table  8 Vibration  Isolator  Design,  Case  (i) 
(Second  Order  Algorithm) 


iter 

*>i 

obj  . 

1 1 5b1 1 1 

1 1 So" ! Imax 

tight 

constraint (0) 

0 

0.15 

1.0 

53.*» 

22.35 

(5)0.0017 

( • 

865) 

1 

0.1416 

0.9713 

43.42 

22.34 

(0)0.0002 

gl(‘ 

857  ) 

2 

0.1207 

0.8998 

29.69 

22.30 

(0)0.0075 

(1 

.087) 

3 

0.1182 

0.9389 

33.94 

19.24 

(1)0.004 

gl(  ' 

637) 

h 

0.1385 

0.9187 

24.06 

15.66 

(0)0.03 

837 ) ,g, (1.113) 

5 

0.1621 

0 . 9084 

21.40 

3.61 

(1)0.0001 

gx(  • 

843)^(1.065) 

6 

0.1686 

0.9090 

21.06 

0.015 

(0)0.0004 

gi(' 

848),e,(1.059) 

7 

0.16860 

0.90906 

21.060 

0.0008 

(0)0.5(-6) 

?1(  • 

648) ,g1(l .059) 

Table  9 Vibration  Isolator  Design,  Case  (ii) 


(Second  Order  Algorithm , y fixed ) 


iter 

bi 

obj  . 

lUbhi 

11  <5 a1  Umax  tight  constraint ,g< (a) 

0 

0.15 

1.0 

37.8 

19.99 

(4)0.0015 

Sx( -9) 

1 

JL. 

0.15 

1.0 

37.8 

19.99 

( 0 ) 0 . 2 (-4  ) 

gx(.9) 

2 

0.15 

1.0 

38.24 

19.99 

(0 ) 0 . 9 (-6 ) 

g ( .9) 

3 

0.15 

1.0 

38.24 

19.99 

(0 ) 0 . 7(-6 ) 

( • 9 ) 

n 

0.15 

1.0 

38.25 

19.99 

(0)0. 5 ( —6 ) 

gx(.9) 

5 

0.15 

1.0 

38.25 

19.99 

(0)0. 5 ( — 6 ) 

( . 9 ) 

6 

0.15 

1.0 

38.25 

0.90 

(0)0 . 5(-6 ) 

g1(.9),g1(1.13) 

7 

0.1575 

0.9537 

22.84 

9.75 

(1)0.001 

-9)  ,g1(l.  092) 

8 

0.1380 

0.9479 

17.98 

6.12 

(1)0.002 

g1(-9),g1(1.108) 

0 

0.1299 

0.9485 

16.65 

3.53 

(0)0 . 6(-4 ) 

g1(.9),g1(1.112) 

10 

0.1263 

0.9516 

16.61 

1.30 

(0)0. 5(-4) 

g1(.9),g1(l.H6) 

11 

0.1250 

0.9527 

16.61 

0.49 

(0)0.6(-5) 

g1(.9),g1(l.H8) 

12 

0.1246 

0.9531 

16.61 

0.18 

(0)0. 3 (-6) 

81( *9) ,g1(l .118) 

13 

0.1244 

0.9532 

16.61 

0.064 

(0)0.2(-6) 

g1(.9),g1(1.119) 

14 

0.12432 

0.95329 

16.605 

0.008 

(0)0.5(-6)g 

g,  ( . 9)  ,gi  (1.119) 

Table  10  ^ive  Degree-of-Freedom  Vehicle  Problem 


(Second  Order  Algorithm) 


iter 

b< 

Ji 

ob  j . 

lUbhi 

maximizing  a's  1 1 6aA | | max 

violated 

inner 

0 

0.2 

0.2 

0.2 

16.741 

1.0 

0.5  1.2 

3.5  7.0 

Si 

1 

0.288 

0.309 

0 . 21 6 

8.259 

1.0 

1.01  2.0 

5.5  (4)0.15 

gl 

2 

0.394 

0.401 

0.236 

5-191 

0.99 

0.99  1.3 

4.8  (4)0.15 

S1  ’ S2 

3 

0.527 

0.502 

0.268 

3.554 

0.99 

0.98  4.0 

(4)0.15 

*1 

4 

0.688 

0.619 

0.315 

2.637 

0.95 

0.96  3.3 

(4)0.15 

51 

5 

0.873 

0.753 

0.385 

2.107 

0.85 

0.96  2.5 

(4)0.15 

Si 

6 

1.067 

0.893 

0.475 

1.804 

0.66 

0.97  1.8 

(4)0.15 

> &2 

7 

1.250 

1.023 

0.572 

1.636 

0.47 

1.03  1.08 

(4)0.05 

gl»  g2 

8 

1.512 

1.185 

0.749 

1.516 

0.11 

1.429 

(3)0.2(-3) 

S1 

9 

1.742 

0.973 

0.885 

1.4  96 

0.19 

1.470 

(0)0.12 

gl 

10 

1.968 

0.938 

1.027 

1.491 

0.06 

1.435 

(1)0.001 

gl 

11 

2.048 

0.926 

1.104 

1.490 

0.004 

1.505 

(0)0.004 

gl 

12 

2.097 

0.913 

1.141 

1.490 

0.004 

1.505 

(1 ) 0 .6 (-5) 

gl 

13 

2.097 

0.914 

1.141 

1.490 

0 . 5 ( - 3 ) 

1.058 

(0)0.003 

gl 

14 

2.0985 

0.9133 

1.1448 

1.4898 

0 . 2 (-4 ) 

1.5053 

(0)0. l(-4) 

gl 

VO 

VO 


Table  11  Chebyshev  Approximation  Problem 


(Second  Order  Algorithm) 


iter 

b1 

ob  .1  . 

jW.II 

1 1 lS.a^  1 1 max 

violated  Ei(a) 

0 

0.7 

0.3 

-0.4 

-2.4 

0.13(-4) 

0 . 6 ( -2 ) 

(5)0.02 

g^l . o) 

1 

0.7 

0.3 

-0.3928 

-2.393 

0.74C-6) 

0.3C-12) 

(3)0.02 

g1(0),g1(0.42) 

2 

0.700 

0.300 

-0.3927 

-2.393 

0.7K-6) 

0 . 9 ( — 2 3 ) 

(5)0.02 

gn  (0)  »Ei  (0.  **17) 
gjd.o)1 

3 

0.700 

0.300 

-0.3928 

-2.393 

0 . 72 (-6 ) 

0 . 3 ( — 36 ) 

(2)0.7 (-4) 

g1(0),g1(0.4l7) 

gl(1.0) 

100 
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. 5 Conclusions 

In  the  present  chapter,  the  sensitivity  analysis  that 
gives  a first  order  feedback  law  and  the  method  of  solving 
the  approximate  problem,  obtained  from  the  PCD  problem, 
have  been  described. 

Several  examples  are  solved  and  compared  with  the 
solutions  obtained  by  the  first  order  algorithm  of  Chapter 
3.  The  second  order  method  exhibits  rapid  convergence, 
once  the  iterates  are  near  the  solution.  Since  the  second 
order  method  accounts  for  the  change  in  the  maximum  point, 
due  to  small  changes  in  the  design  and  compensates  any 
error  in  the  current  maximum  point,  the  method  does  not 
require  separate  solution  of  the  inner  problems. 

The  disadvantage  of  the  method  lies  in  the  determina- 
tion of  a good  initial  estimate  and  solution  of  the  Kuhn- 
Tucker  equations.  The  method  is  not  reliable  unless  the 
problems  under  consideration  have  nice  properties  such  as 
convexity,  at  least  locally,  or  unless  the  initial  estimates 
are  close  enough  to  the  solution. 

To  make  use  of  the  advantages  of  the  first  and  second 
order  algorithms,  some  variations  will  be  considered  in  the 
following  chapter. 
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CHAPTER  5 
HYBRID  METHODS 

5.1  Introduction 

In  the  previous  chapters,  purely  first  and  second  or- 
der algorithms  were  discussed.  It  has  been  observed  that 
one  numerical  approach  is  not  always  effective  for  all 
problems.  There  may  be  a variation  on  these  approaches 
that  is  better  suited  to  a given  problem.  It  is  the  pur- 
pose of  this  chapter  to  discuss  some  possible  variations 
of  the  algorithms  described  in  the  preceding  chapters. 

Once  one  defines  the  effectiveness  of  an  algorithm, 
in  terms  of  rate  of  convergence  and  the  effort  of  con- 
structing, and  running  a program,  a best  combination  may 
be  possible  for  a given  problem.  In  this  sense,  the  se- 
lection of  approach  is  an  optimization  process,  which  is 
critical  in  practical  problems.  In  the  following,  possible 
alternatives  are  given  for  the  POD  algorithm  and  a com- 
parison of  computational  effort  is  shown,  based  on  the 
number  of  numerical  operations  required  to  implement  the 
algorithm. 
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5.-JL  Review  of  Alternatives 

As  stated  earlier,  the  first  order  algorithm  was 
developed  under  the  assumption  that  the  subproblem  is 
solved  as  an  integral  part  of  each  iteration.  Hence,  one 
has  freedom  to  choose  any  solution  method  in  solving  the 
subproblem.  In  the  second  order  algorithm,  the  solution 
of  subproblem  is  adjusted  through  use  of  sensitivity  fac- 
tors, simultaneously  with  the  minimization  of  the  objective 
function.  Also,  a first  order  feedback  law  was  developed, 
with  a correction  term  for  any  error  in  the  solution  of  the 
subproblem.  A full  second  order  approximation  in  the  de- 
sign variable  was  necessary  to  maintain  the  order  of  the 
approximation.  From  a numerical  point  of  view,  however,  a 
first  order  method  such  as  steepest  descent  has  been  found 
to  be  effective  for  the  solution  of  NLP.  It  is  natural  to 
attempt  a first  order  outer  procedure,  with  sensitivity 
information  obtained  from  a second  order  approximation  of 
the  Inner  problem. 

Variations  in  methods  fall  into  two  categories:  One 

is  essentially  a first  order  algorithm  with  a higher  order 
method  of  solution  of  the  subproblem.  The  other  is  a vari- 
ation from  the  second  order  algorithm  and  uses  a first  or- 
der approximation  In  the  outer  procedure.  Since  this 
method,  then,  uses  the  first  order  total  derivatives  with 
respect  to  b as  noted  in  Section  4.3,  this  method  is. 
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in  reality,  first  order.  Including  the  original  first  and 
second  order  algorithms,  four  alternative  methods  are 
possible.  They  are: 

(1)  first  order  maximization  + first  order  outer  procedure, 

(2)  second  order  maximization  + first  order  outer  proce- 
dure , 

(3)  sensitivity  with  second  order  maximization  + first 
order  outer  procedure,  and 

(4)  sensitivity  with  second  order  maximization  + second 
order  outer  procedure. 

The  two  alternative  methods,  (1)  and  (2),  have  their 
background  in  the  analysis  given  in  Chapter  3,  where  it  is 
shown  that  the  first  order  algorithm  for  the  solution  of 
the  POD  problem  is  essentially  an  alternation  of  maximiza- 
tion and  minimization,  as  is  frequently  used  in  min-max 
problems.  The  alternative  methods  (3)  and  (4)  do  not 
require  a separate  procedure  for  soving  inner  problems, 
since  the  sensitivity  analysis  with  error  compensation, 
presented  in  Chapter  4 , gives  the  necessary  change  in  the 
solution  point.  Being  first  order  in  nature,  method  (3) 
is  expected  to  have  sufficient  reliability,  while  maintain- 
ing the  accuracy  of  maximum  points.  On  the  other  hand,  the 
full  second  order  method  (4)  is  too  complicated  to  solve 
the  resulting  approximate  problem,  requiring  a Newton's 
procedure . 
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5.3  Second  Order  Method  of  Maximization 
If  one  puts  6b=0  in  the  sensitivity  analysis  of  Chap- 
ter one  has  a second  order  method  of  maximization.  In 
this  section,  this  approach  is  examined  more  closely,  since 
it  can  be  implemented  for  NLP  and  can  be  used  in  method 
(2). 

Since  it  has  been  assumed  that  the  a-conctraints  are 
linear  in  a,  the  problem  considered  here  is  to  maximize 
AT6a  + ^6cT  AA«a  (5.1) 

subject  to 

Q^6a  < Aq,  (5 .2) 


and 


T ? 

6a  Wa6a  £ £ . 

(5.3) 

From 

Section  2 . ^ , the  solution  of  the 

Kuhn-Tucker  necessary 

conditions  is  obtained  as 

1 2 

6a  = -6  a + 6a  y 

(5.4)- 

where 

6a1  = ♦ 1 ( I-M  ) A , 

(5.5) 

6a?  = ♦ -1M  Aq 

q 

(5.6) 

and 

* = AA  - 2vW  , 
a 

(5.7) 

M = Q(QT*  _1Q  )-1 , 

(5.8) 

T -1 

M = M Q ♦ . 

qq  q 

(5.9) 

As  direct  multiplication  shows,  M is  a projection  ma- 

QQ 

trix : i.e.,  M M =M  . To  interpret  this  as  in  the 
qq  qq  qq 


106 


gradient  projection  method,  the  following  relations  can  be 
verified : 


T 


(1) 

1 2 

6o  *6a  * 0, 

(5.10) 

(2) 

T 2 

Q 6a  = Aq , 

(5.1D 

(3) 

T 1 

Q 6ax  = 0, 

(5.12) 

(4) 

T 1 1 lT  1 

A 6a  + p-6 a AA6a  0, 

(5.13) 

where  $ is  negative  semi-def inite . 

Property  ( 4 ) , combined  with  the  yet  undetermined  con- 
stant is  valuable  in  practical  computation.  By  choosing 
v>0  large  enough  such  that  4 is  negative  definite,  5a^  is 
always  in  the  direction  of  increasing  the  objective  func- 
tion (5.1).  Hence,  one  can  generate  a sequence  of  points 
converging  to  a local  maximum  point.  This  may  lead,  how- 
ever, to  a false  impression  of  convergence,  if  one  chooses 
very  large  v . 

Much  literature  has  been  published  on  Newton's  method 
and  its  variations  [45,46],  Most  of  this  literature  treats 
unconstrained  problems.  When  constraints  are  absent,  the 
method  presented  here  is  another  version  of  Newton's  meth- 
od. Same  techniques  have  been  studied  in  the  literature 
[45,47],  from  a different  point  of  view.  Although  the  re- 
liability of  the  method  can  be  easily  guaranteed,  conver- 
gence properties  are  largely  dependent  on  the  optimal  choice 
of  v [453.  An  estimate  of  v and  numerical  implementation 
are  discussed  in  Section  6.3. 
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5 . 4 A Computational  Modification 
for  Alternative  Method  (3) 

As  has  been  discussed  in  Section  5.2,  method  (3) 
is  rather  efficient.  In  this  section,  a computationally 
efficient  modification  is  described  for  the  case  in  which 
no  a-constraints  are  present,  i.e.,  M ®0=M  in  Eq.  (4.7). 

'i  MM 

In  the  case  of  the  first  order  outer  procedure,  the 
inversion  of  ♦ for  the  expression  (4.14)  can  be  avoided  in 
the  following  way.  From  Eq.(4.7),  multiplying  by  ♦, 

♦ 6a  = -D'6  b - C'  (5.14) 

where , 

♦ = AA  - 2vW  , 

a * 

T 

D'-  BA  , 

C'=  A. 

Method  (3)  utilized  a first  order  expansion  in  Eq . (4.4), 
6g  = BT6b  + AT6a . (5.15) 

To  eliminate  the  second  term,  let  < be  determined  as  the 
solution  of 

♦ T«  = A.  (5.16) 

Post-multiplying  by  6a,  after  taking  the  transpose,  one 
obtains 

tc  T^6  a = AT6  a . (5.17) 

Hence , 

6 g = (BT-<TD')6b  - <V.  (5.18) 

Therefore,  the  first  order  approximate  problem,  using 
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sensitivity  Information,  Is  obtained  as:  Minimize 


6 f 

S 

tT 

iJ  fib 

(5.19) 

subject 

to 

. T 

ae 

i * fib  - K+  <_  Ag 

(5.20) 

where 

3 f T 

(5.21) 

3b  * 

£J 

= 

(B  - D'TOj, 

(5.22) 

ss 

-(<V)j,  J e Ic. 

(5.23) 

Here,  the  subscript  j denotes  the  quantity  calculated  for 
the  J-th  inner  problem.  After  solving  this  approximate 
problem  for  fib,  using  the  gradient  projection  method  de- 
scribed In  Section  2.4,  6a  is  obtained  from  Eq.  (5.1*0. 

It  is  noted  that  the  inversion  of  ♦ has  been  replaced 
by  two  solutions  of  a system  of  linear  equations  (5-16) 
and  (5.14),  which  i3  a great  computational  advantage  when 
the  dimension  of  a is  large. 

5.5  An  Application  to  the  Optimal  Design 
Problem  with  State  Equations  Given  Implicitly 
by  a Minimum  Principle 

The  type  of  problems  under  consideration  has  been  in- 
troduced in  Section  2.2.3,  and  can  be  treated  by  the  meth- 
od discussed  in  Section  5.4.  In  the  present  section,  the 
case  of  a function  H(b,z)  in  Eq.  (2.15),  which  is  quadratic 
with  respect  to  z , is  considered.  Comparisons  are 
possible,  since  this  results  in  a common  optimal  design 
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problem.  Most  of  the  problems  in  mechanics  belong  to  this 
category,  where  H(b,z)  usually  represents  total  potential 
energy  of  the  system.  The  problem  is  stated  as:  Minimize 

f(b,z)  (5.24) 


subject  to 

g(b  ,z)  <_  0,  (5.25) 

where  z maximizes  the  quadratic  function, 

-H ( b , z ) = ^(2zTF  - zTK(b)z),  (5-26) 

K(b)  is  a positive  definite  matrix,  and  F is  a known 
vector. 

The  sensitivity  analysis  of  Section  4.2,  with  the 
modification  given  in  the  previous  section,  is  applied  to 
the  present  problem.  Noting  that  the  state  variable  z 
corresponds  to  the  free  parameter  o,  one  has 

♦iz  = -D'fib  - C'  (5.27) 


where , 


* * -K  - 2vW, 

D'«  - |^[K(b)z], 

C'=  F - Kz. 

With  these  interpretations,  the  optimal  design 

comes,  from  the  linearized  expressions  of  Eqs. 

(5.25),  minimize 
JT 

6 f 6b 


(5.28) 

(5.29) 

(5.30) 
problem  be- 
(5.24)  and 

(5.31) 


subject  to 


♦ T 

l 6b 


K*  i &g. 


(5.32) 
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whe  re 


(5.33) 

= 3KT  _ d't** 

3b 

<5.3*i  ) 

K*  = -k*c'  , 

and  **  are  solutions  of  the  systems 

(5.35) 

.T  J _ L£t 
* * " 3Z  > 

(5.36) 

♦V  = 

3 z 

(5.37) 

This  linearized  approximate  problem  can  be  solved  by  the 
gradient  projection  method  described  in  Section  2D.  Then, 
the  state  variable  z is  predicted  as  z+jz,  where  4 z is 
solved  from  Eq.  (5-27). 

Note  that  the  state  equations  implicitly  given  by  Eq . 
(5. 2C)  can  be  obtained  explicitly  as 

K(b)z  - F = 0.  (5.38) 

The  final  linearized  approximate  problem  of  Eqs . (5. 3D 
and  (5.32)  Is  essentially  the  same  as  the  one  obtained  di- 
rectly using  Eq . (5-38)  [33].  In  the  present  approach, 

4> 

the  problem  contains  new  terms  -2vW  and  K , which  give  an 

additional  feature  to  the  method.  The  factor  v,  acting  as 

damping,  stabilizes  the  convergence  of  the  method,  as  ex- 

4> 

plained  in  Section  5.3,  and  K compensates  the  error  in- 
volved in  the  current  solution  z. 

In  the  case  of  a quadratic  function  H(b,z),  little 
computational  efficiency  is  gained  but  in  the  general  case 
of  state  relations  given  implicitly  by  a nonlinear  pro- 
gramming problem,  no  other  method  seems  known. 
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5.6  Numerical  Examples 

The  weapons  allocation  problem,  solved  in  previous 
chanters  by  the  first  and  second  order  methods,  is  solved 
by  the  hybrid  method  (3).  Table  12  shows  that  the  hybrid 
method  has  the  advantage  of  the  second  order  algorithm, 

In  automatic  adjustment  of  maximum  noints,  and  it  has  the 
reliability  of  the  first  order  algorithm.  Because  the  max- 
imum points  are  self-adjusted,  better  convergence  is  obtain 
ed  than  the  pure  first  order  algorithm.  The  computer  time 
for  this  problem  was  about  0.1  seconds  per  Iteration. 

Application  of  the  hybrid  method  (3)  to  the  vibration 
isolator  design  also  shows  that  the  method  can  be  superior 
to  the  other  two.  The  computer  time  per  iteration  was 
0 . 1 seconds . 

The  solution  of  the  5 degree-of-freedom  vehicle 
problem  is  shown  in  Table  15-  Although  at  the  beginning 
the  convergence  was  rapid,  it  was  slow  near  the  solution. 
The  computer  time  in  this  case  was  0.7  seconds. 

Results  for  the  Chebyshev  problem  in  Table  16  show  a 
similar  trend  of  rapid  convergence  for  the  first  few  itera- 
tions. Although  the  mathematical  solution  is  difficult  to 
attain  by  the  hybrid  method,  since  this  method  is  essential 
ly  first  order,  the  solution  obtained  seems  good  for  engi- 
neering purposes.  The  computer  time  per  iteration  was  0.23 
seconds,  and  Is  about  'J0%  of  that  for  the  second  order 
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Table  12  Weanons  Allocation  Droblem 
(Hybrid  Method) 

, tieht 

iter  ob.1  . I Ub1 1 | ^i IUa  j|  Si,  Qj 

0.1  0.2  0.3 


0 

0.1 

0.2 

0.3 

-2.10 

0.795 

0.0 

C.2H6 

0.563 

(4)0.002 

ql 

1 

0.09  9 

0.200 

0.300 

-2.10 

0.792 

0.0 

0.21414 

0.567 

(0)0.0006 

ql 

2 

-0.0014 

0.145 

0.245 

-2.22 

0.465 

0.0 

0.127 

0.  Hi)  14 

(0)0.02 

3 

0.0 

0.327 

0.351 

-2.03 

0.870 

0.0 

0.305 

0.536 

(4)0.007 

n 

14 

0.0 

5.19^ 

0.245 

-2.230 

0.180 

0.0 

0.19*1 

0. 14115 

(5)0.004 

n 

5 

0.0 

0.181 

0.288 

-2.235 

0.077 

0.0 

0.183 

0 . *15*1 

(0)0.014 

,1  ! 

6 

0.0 

0.172 

0.281 

-2.23*1 

0.062 

0.0 

0.167 

0.H71 

(3)0.002 

11 

7 

0.0 

0.179 

0.279 

-2.23*1 

0.03*1 

0.0 

0.169 

0.1171 

(0)0.002 

n 

8 

0.0 

0.182 

0.274 

-2.23*1 

0.016 

0.0 

0.172 

0.470 

(0)0.001 

1 

q 

0.0 

0.183 

0.273 

-2.23*1 

0.01 

0.0 

0.173 

0.469 

(0)0.001 

n 

10 

0.0 

0.183 

0.272 

-2.23*1 

0.006 

0.0 

0.17*1 

0.468 

(0)0.001 

it 

11 

0.0 

0.18277 

0.27174 

-2.2340 

0.0025 

0.0 

0.17*19 

0.4663 

(1)0.0003 

n 

£TT 


Table  13  Vibration  Isolator  Design,  Case  (1) 
(Hvbrid  Method) 


iter 

bi 

ob  j . 

1 la  b1 1 1 

I I ! Imax 

tight  constraint ,g(a) 

0 

0.15 

1.0 

58.4 

22.35 

(5)0.0017 

g1( . 865) 

1 

0.1507 

0.9447 

30.84 

22.31 

(0)0.006 

OO 

-tr 

VO 

2 

0.1511 

0.9057 

22.78 

9.04 

(1)0.0018 

( • 835) ,E1 ( 1 . 074  ) 

3 

0.1636 

0.9084 

21.18 

3.44 

(0)0.01 

g1(.835),g1(1.074) 

4 

0.1761 

0.9082 

21.19 

4.38 

(1)0.0001 

g1(.851),g1(1.054) 

5 

0.1636 

0.9095 

21.12 

2.88 

(0)0.0016 

gx(.  843)^(1.065) 

6 

0.1685 

0.9091 

21.07 

0.07 

(0)0.0004 

gx(.848) ,g1(1.059) 

7 

0.1686 

0.9091 

21.06 

0.03 

(0)0.2(-5) 

gl(.848),gl(1.059) 

8 

0.16857 

0.90906 

21.060 

0.02 

(0)0. 8 ( — 6 ) 

g1(.848),g1(1.059) 

'lable  1 4 Vibration  Isolator  Design,  Case  (11; 
(Hybrid  Method) 


iter 

bi 

ob.l  . 

| |5b 

0 

0.15 

1.0 

37.8 

19. 

1 

0.1494 

0 . 9t55 

24.27 

is. 

2 

0.11487 

0.9443 

19.25 

19. 

3 

0.1480 

0.9309 

17.89 

11. 

14 

0.1310 

0.9472 

16.75 

4. 

5 

0.1161 

0.9602 

16.73 

5. 

6 

0.1310 

0.9475 

16.69 

4. 

7 

0.1223 

0.9550 

16.62 

1. 

8 

0.1271 

0.9509 

16.62 

1. 

9 

0.1251 

0.9527 

16.61 

0. 

10 

0.1247 

0.9530 

16.61 

0. 

11 

0.12414 

0.9532 

16.61 

0. 

12 

0.1214 14 

0.9533 

16.61 

0. 

13 

0.12432 

C. 95329 

16.605 

0. 

I 1 6 a " . 1 max 

tight  constraint ,£i ) 

(4)C. 0015 

P]/  .9) 

(0)0.0015 

g1( . 9) 

(0)0.0003 

^(•9) 

(0)0.0002 

Px ( . 9 ) ,g1(1.086) 

(0)0.0015 

g1(.9),g1(l.H2) 

(0)0.0007 

K1(.9),g1(1.13),g2(.94) 

(0)0.001 

g1(.9),g1(1.112) 

(0)0.0002 

g1(.9),g1(1.12),g2(.94) 

(0)0.0001 

g1(  .9)  jg^d.ns) 

(0)0. l(-4) 

P1(. 9), gx(l. 118) 

(0)0.1  (-6) 

g1(.9),g1 (1.116) 

(0)0.4 (-6) 

g1(.9),g1(l.H8) 

(0)0. 7(-6) 

g1(.9),g1(1.119) 

(0 ) 0 . 6 ( -6 ) 

g1(.9),g1(1.119) 

ill 

99 

95 

89 

53 

28 

63 

20 

3*< 

85 

50 

27 

09 

03 

005 
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Five  Derree-of-Freedor.  Vehicle  Probler. 


(Hybrid  Method) 


ter 

h 

obj. 

lUbMl 

! 1 <5 a-1- j | max 

violated  g-»  (a ) 

0 

0.2 

0.2 

0.2 

16.741 

1.0 

(4)0.15 

p-j^  (0.5) 

1 

0.288 

0.309 

0.216 

8.259 

1 . 0 

(4)0.15 

p1(1.01) 

2 

0.39^ 

0.401 

0.236 

5.191 

0.99 

(4)0.15 

P1 ( 0 . 95 ) ,£1 (1 . 3) 

3 

0.527 

0.502 

0.268 

3.554 

0.99 

(4)0.15 

P1 ( 0 . 9 S ) 

4 

0.688 

0.619 

0.315 

2.637 

0.95 

(4)0.15 

g1(0.96) 

5 

0.873 

0.753 

0.385 

2.107 

0.85 

(4)0.15 

P1 C 0 . 9 6 ) 

6 

1.068 

0.8$4 

0.475 

1.803 

0 . 66 

(4)  0.002 

g1(0.97) 

7 

1.251 

1.024 

0.572 

1.635 

0.47 

(1)0.003 

^(1.02?) 

8 

1 . 396 

1.124 

~ 

0 • s l 

1.551 

0.31 

(0)0.047 

P1(1.09S) 

9 

1.493 

1.1&5 

0.72  0 

1.520 

0.15 

(2)0. 3(-3) 

( 1 . 236 ) 

10 

1.536 

1.195 

0.755 

1.514 

0.096 

(0)0.06 

e1(i.496) 

15 

1.067 

1.110 

0.638 

1.504 

0.062 

(0)0.002 

g1(1.484) 

20 

1.740 

1.052 

0.942 

1.497 

0.14 

(0)0.036 

g1(1.593) 

25 

1.775 

0.657 

0.974 

1.466 

0.08 

(5)0.04 

(1. 633) 

30 

1.838 

0. 9S9 

1.006 

1.492 

0 . 031 

(0)0.016 

gx  (1 . 493) 
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Table  It*  Chebyshev  Approximation  Problem 
(Hybrid  f'ethod) 


iter 

b, 

A. 

obj  . 

II^IL 

1 I 6a1 1 I max 

violated  gi(a) 

C 

0.7 

0.3 

-0.4 

-2.4 

0 . 1 3 ( -4  ) 

0 . 6 ( -2 ) 

(5)0.02 

g;L(1.0) 

1 

0.7 

0.3 

-0.4 

-2.4 

0.13C-4) 

0 . 6 (-2 ) 

( 5) 0 . 5(-4 ) 

gx(1.0) 

2 

0.7326 

0.3044 

-0.377 

-2.399 

0 . 14 ( -2  ) 

0 . 3 ( -2 ) 

(0)0.2 (-8) 

gl(o.),gl(i.o) 

3 

0.7038 

0.2993 

-0.4046 

-2.399 

0 . 10(-4  ) 

0 . 3 ( -2 ) 

(0)0. l(-8) 

gl(0. ) ,gl(l .0) 

h 

0.7060 

0.2930 

-0.3978 

-2.399 

0.1K-5) 

0.8(-3) 

(0)0 ,2(  — £) 

g1 ( 0 . ) ,gl(l. 0) 

5 

0.7058 

0.2937 

-0.3988 

-2.399 

0.27 (-6) 

0 . 6(-3) 

(0)0.0 

E1(0.),g1(1.0) 

M 
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5.7  Comparison  of  Computational  Effort 
In  this  section,  the  effect  of  the  dimension  on  the 
overall  number  of  arithmatic  operations  will  be  given,  in 
an  approximate  way.  The  count  given  denotes  the  number  of 
mulltplications . Nearly  the  same  number  of  operations  of 
addition  or  subtraction  Is  necessary.  Operations  such  as 
division  and  logic  are  not  included  here. 

The  computational  effort  of  one  iteration,  for  a full 
second  order  algorithm,  is  proportional  to 

s with  a factor  of  (n^m^  + nm2p  + n^p^* • ♦ ) , 

plus  (STATE,  G,  MATINV(p) , GAUSS (p ) , Q,  MATINV(g) } , 
nc  with  a factor  of  (sm  + 4g» • • } , 

with  a factor  of  (sCn^  + np  + p^)***), 
p^  with  a factor  of  (s(n  + q)***}, 
n‘  with  a factor  of  (s(n^  + m^  + q^)***), 
and  g^  with  a factor  of  (s  + 2n***}.  Also,  GAUSS(n+g)  is 
necessary  as  many  times  as  the  iterations  of  Newton's 
method . 

In  the  above,  n,  m,  and  p denote  the  number  of  design 
variables,  state  variables,  and  free  parameters,  respec- 
tively. The  numbers  s,  g,  and  q are  the  number  of  inner 
nroblems,  active  design  variable  constraints,  and  active 
free  parameter  constraints,  respectively;  STATE  solves  the 
state  equations  and  calculates  necessary  derivatives;  G 
and  0 denote  the  function  evaluations  of  g's  and  q'r, ; 
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NIATINV(p)  denotes  a p-dimensional  matrix  inversion;  and 
GAUSS(p)  denotes  a p-dimensional  linear  equation  solver. 

The  notation  Implies  lower  order  terms  are  neglected, 
where  n > m > p is  assumed. 

In  the  case  of  the  first  order  algorithm,  the  compu- 
tational effort  of  one  iteration  is  proportional  to  s,  with 
a factor  of  (nm  + (m+2n)p  + g^  • • • } plus  {STATE,  G,  Q, 
GAUSS(g)  },  and  r2  with  a factor  of  {s  + 2n  • • • } . 

It  is  noted  that  the  Increase  of  computational  effort 
for  the  second  order  algorithm,  compared  with  that  of  first 
order,  is  enormous  and  it  depends  on  the  convergence  of 
the  Newton's  method.  For  the  hybrid  method  (3)>  the  only 
increase  in  the  computational  effort  will  appear  in  the 
solutions  of  inner  problems. 

5.8  Conclusions 

In  this  chapter,  modifications  of  the  algorithms 
described  in  previous  chapters  have  been  presented.  The 
hybrid  method  (3),  which  is  essentially  of  first  order,  has 
been  found  quite  effective.  This  is  because  it  makes  use 
of  both  the  advantage  of  sensitivity  information  and  the 
reliability  of  first  order  outer  procedure. 

A second  order  method  of  solution  for  a nonlinear 
programming  problem  has  been  discussed  in  detail,  which  has 
appeared  in  the  sensitivity  analysis  in  Section  *J.2.  Also 
an  application  of  sensitivity  analysis  to  ar.  optimal  design 
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problem,  where  state  equations  are  plven  implicitly  by  a 
minimum  principle,  is  discussed. 
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CHAPTER  6 

IMPLEMENTATION  AND  PROGRAMMING 
6_.  1 Introduction 

In  this  chapter,  the  description  of  a computer  pro- 
gram, based  on  the  second  order  algorithm  developed  In 
Chapter  , Is  given.  The  program  for  the  first  order  algo- 
rithm can  be  obtained  by  removing  all  the  statements  In- 
volving second  order  derivatives.  A schematic  flow  chart 
of  the  second  order  algorithm  follows  a general  description 
of  the  program.  In  the  following,  Newton's  method  may  be 
replaced  by  the  gradient  projection  method  for  the  first 
order  outer  procedure,  i.e.,  methods  (1),  (2),  and  (3). 

The  other  modifications  are  self-evident  from  Chapter  5* 

The  algorithms  discussed  in  the  previous  chapters  are 
rather  complicated.  This  is  due  partly  to  the  complexity 
of  the  problem  formulation,  but  also  to  the  complexity  of 
the  underlying  method,  i.e.,  the  gradient  projection  method 
which  needs  delicate  computational  skills. 

In  this  chapter,  the  second  order  algorithm  Is  conven- 
iently divided  into  several  subalgorithms,  or  subprocedures. 
Unfortunately,  most  of  subprocedures  work  nicely  only  for 
problems  that  have  special  properties. 
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In  the  present  thesis,  the  algorithm  is  motivated  by 
engineering  purpose  and  is  numerical  solution  oriented. 
Hence,  the  major  contents  of  the  present  chapter  are  based 
on  numerical  experience  with  a number  of  engineering 
problems  and  contains  the  description  of  parameter  choice 
for  implementation,  a discussion  of  difficulties  Involved, 
and  an  inclusion  of  a heuristic  step  for  certain  types  of 
problems.  As  is  true  with  the  rather  well  established 
methods  of  solving  nonlinear  programming  problems,  such  as 
the  gradient  projection  method,  method  of  feasible  direc- 
tions, and  various  other  steepest  descent  techniques, 
there  is  a trade-off  between  being  mathematically  precise 
and  reducing  the  amount  of  computational  and  programming 
effort . 


6.2  Description  of  Computer  Program 
6.2.1  Simplified  Block  Diagram 
of  POD  Algorithm 

As  shown  in  Figure  7,  the  whole  program  consists  of 
five  major  steps.  With  the  estimated  designs  and  para- 
meters, derivatives  and  function  values  are  calculated  by 
subroutines.  After  checking  the  constraints,  various 
matrices  involved  in  the  theoretical  algorithm  of  the  pre- 
vious chapter  are  constructed.  If  necessary,  a maximiza- 
tion procedure  Is  executed,  using  the  matrices  calculated 
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Figure  7 Simplified  Block  Diagram  of  POD 
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above.  The  Newton's  method  is  applied  to  solve  for  design 
modifications  and  Lagrange  multipliers.  As  is  the  case 
in  the  algorithm  developed  in  the  previous  chapter,  signs 
of  the  multipliers  corresponding  to  inequality  constraints 
must  be  checked.  If  some  are  negative,  associated  con- 
straints are  eliminated  and  the  active  constraint  set  is 
redefined.  A stopping  criteria  is  checked  and  iteration 
is  terminated,  or  is  continued  with  improved  designs  and 
free  parameters.  In  the  following  sections,  these  major 
steps  are  described  in  detail. 

6.2.2  Explanation  of  the  Major  Steps 
(1)  If  state  equations  are  involved,  a subroutine 
(STATE)  is  constructed  to  solve  for  the  values  of  state 
variables  and  the  matrices  described  in  Appendix  D.  Sub- 
routines are  constructed  for  the  functionals  of  subprob- 
lems, in  order  to  construct  the  necessary  derivatives  and 
functional  values.  At  this  stage,  the  constraint  function- 
als are  checked  for  violations.  Even  though  a parametric 
constraint  is  strictly  satisfied,  it  is  necessary  to  re- 
tain all  derivatives  required  to  adjust  the  worst  case, 
since  the  maximum  point  is  going  to  be  changed  in  later 
steps,  and  may  become  an  active  constraint.  It  is  observed 
in  nonlinear  programming  that  if  the  iterates  are  near 
the  solution,  the  active  set  remains  the  same  during  the 
iterations.  In  this  case  the  inactive  set  can  be  dropped. 


eliminating  unnecessary  computations.  Caution  is  necessary 
to  retain  the  global  maximum  of  the  subproblem,  when  it  is 
not  a concave  programming  problem. 


(2)  With  the  return  of  derivatives  and  function  val- 
ues from  the  appropriate  subroutines,  the  various  matrices 
in  Eqs.  (4.15),  (4.16),  (4.17),  and  (4.7)  are  constructed. 

In  constructing  the  matrices  and  L1*1  in  Eq.  (4.14),  the 

active  constraints  are  defined  through  use  of  an  e-active 

Index  set,  Introduced  in  Section  2.4.2.  By  removing  all 

T l 

the  inequalities,  g^  + 6b+  ^b^L^fib-K^  1 0 > for  which 

< 0,  a zigzagging  phenomenon  would  take  place.  To  re- 
duce this  phenomenon,  constraints  with  sufficiently  nega- 
tive values  are  removed,  keeping  weakly  negative  and  strict- 
ly violated  constraints  in  the  active  set.  The  same  idea 
can  be  used  for  a-constraints . 

With  these  index  sets  defined,  the  derivative  matrix 
of  Eq . (4.23)  is  constructed. 

(3)  In  solving  the  algebraic  Eqs.  (4.20),  (4.21),  and 

(4.22)  for  6b^ , and  y by  the  Newton's  method,  an  ini- 

tial estimate  is  made  according  to  the  suggestion  in  Section 

4.3.  Prom  Eq.  (4.24),  | 1 6b1 1 | = 6blTW  6b1  is  calculated, 

w K U 


and  used  as  a stopping  criterion.  Because  of  the  nonline- 
arity of  the  equations,  no  unique  solution  is  expected. 

Among  those  solutions  that  exist,  the  solution  giving  the 

T 

smallest  value  of  6b  W^6b  is  chosen. 
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(4)  With  the  designs  fixed,  a "sufficiently"  accurate 
maximum  point  for  each  subproblem  is  often  desired,  during 
the  iterations  as  well  as  at  the  beginning  of  the  optimi- 
zation process,  even  though  the  second  order  algorithm 
achieves  maximization,  as  Eq . (4.7)  shows.  For  this  pur- 
pose one  may  go  through  the  outer  iterations  to  compute 
the  vector  C of  Eq.  (4.7),  without  calculating  irrelevant 
matrices.  This  step  is  denoted  as  the  maximization  proce- 
dure in  Figure  7.  The  expression  for  the  vector  C in  Eq. 
(4.7)  is  shown  to  be  different  from  that  obtained  from  the 
first  order  gradient  projection  method  (See  Section  2.4), 
only  in  the  matrix  ♦ of  Eq.  (4.11).  This  similarity  arises 
from  the  fact  that  the  a-constraint  was  expanded  only  up 

to  first  order  terms,  in  Eq.  (4.5).  The  algebraic  signs 
of  the  multipliers  associated  with  the  a -constraints , 
given  in  Eq . (4.10),  must  be  checked. 

The  estimation  of  the  undetermined  constant  v,  which 
determines  stepsize,  will  be  considered  in  the  next  sec- 
tion . 

(5)  The  solution  of  the  algebraic  equations  (4.20), 
(4.21),  and  (4.22),  may  not  be  the  desired  solution  of  the 
approximate  problem.  One  has  to  satisfy  the  positivity 
conditions  of  the  multipliers.  Hence,  the  next  step  is  a 
check  of  the  signs  of  these  multipliers,  including  the 
multipliers  corresponding  to  the  a-constraints . 
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Since  these  multipliers  are  implicitly  related,  one  is  not 
sure  whether  the  removal  of  any  a -constraints  correspond- 
ing to  negative  multipliers  can  make  the  negative  multi- 
pliers corresponding  to  the  b-constraints  positive.  At 
this  stage,  the  question  is  not  settled. 

In  the  present  program,  the  multipliers  corresponding 
to  the  a -constraints , u,  are  checked  first.  If  some  are 
negative,  removal  of  the  a-constraints  corresponding  to 
these  negative  multipliers  follows.  At  the  same  time,  the 
sign  of  multiplier,  A,  corresponding  to  b-constraints  is 
checked.  Next,  if  y is  negative,  the  stepsize  constraint 
is  removed. 

(6)  With  all  multipliers  nonnegative,  the  design  and 
free  parameters  are  updated  by  adding  the  increment  formed 
from  the  Newton  procedure,  and  Eq.  (4.7). 

By  the  solution  b and  of  POD,  one  means  that  for 
fixed  b,  is  the  maximum  point  of  g^ , under  the  a-con- 
straints. Hence,  ||fia^||  should  be  zero  (See  Section  2.*J). 
For  fixed  a ^ , b is  the  minimum  point  of  the  objective  func- 
tion subject  to  the  b-constraints,  satisfying  || fibril  =0. 
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6.3  Various  Other  Devices 
6.3.1  The  Choice  of  the  Multiplier  v 
In  Eq.  (4.11),  the  constant  v must  be  known  to  obtain 
the  matrices  D and  C.  It  is  known  that  in  the  first  order 
gradient  projection  method,  the  multiplier  associated  with 
the  stepsize  restriction  is  inversely  related  to  the  step- 
size  in  the  direction  of  steepest  descent  [33].  Like  the 
other  multipliers,  v could  be  determined  in  the  first  or- 
der method  by  using  the  stepsize  restriction,  considered 

as  equality  instead  of  inequality.  The  stepsize  restric- 

2 

tion,  however,  has  an  undetermined  constant  £ , yet  to  be 
decided.  In  practice,  one  chooses  the  stepsize  v and  lets 

p 

C be  determined  accordingly. 

In  the  first  order  gradient  projection  method,  since 

v is  the  inverse  of  the  stepsize,  it  is  sometimes  easy  to 

estimate  a value  for  v based  on  physical  considerations. 

For  the  present  formulation,  however,  since  v does  not 

have  an  explicit  physical  interpretation,  one  must  choose 
2 

S in  Eq.  (4.6)  and  have  v determined  correspondingly. 

The  situation  here  does  not  allow  one  to  compute  v exactly. 
An  approximate  estimation  is  suggested. 

Consider  the  second  order  maximization  described  in 
Section  5.3.  From  Eq . (5.4),  whatever  value  of  v is  chosen 
(assuming  that  this  can  assure  a small  step, 50),  the  term 
compensating  the  constraint  error,  ♦~^Mq&q,  is  determined 
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accordingly . By  this  observation,  let  the  value  of  v be 
determined  with  the  a-constraints  neglected.  Even  in  this 
unconstrained  case,  the  optimality  of  this  multiplier  as 
related  to  stepsize  is  conditional,  as  long  as  the  choice 
of  is  not  optimal.  The  value  of  this  multiplier,  how- 
ever, can  indicate  the  order  of  magnitude  of  the  multiplier 
for  the  subproblem  with  the  a-constraints  included. 

A procedure  for  finding  the  multiplier  for  the  sub- 
problem  with  a-constraints  neglected  follows.  With  6b«0, 
this  procedure  is  exactly  that  of  determining  the  constant 
for  the  Newton' 3 method  of  maximization,  described  in  Sec- 
tion 5.3.  It  was  noted  that  the  approach  there  gives  a 
different  point  of  view  from  those  in  the  literature  C^5» 
31*].  In  this  literature,  v is  determined  so  that  the  com- 
bined matrix  ♦ has  the  property  of  definiteness,  for 
convergence.  Hence,  an  eigenvalue  analysis  or  decomposi- 
tion of  matrices  is  necessary  for  the  determination  of  the 
constant.  In  Section  5-3,  it  has  been  shown  that  v can  be 
interpreted  as  the  Lagrange  multiplier  for  the  stepsize 
restriction.  Hence,  in  the  following,  v will  be  determined 
in  terms  of  s as  a Lagrange  multiplier. 

For  problems  that  are  linear  in  6a,  one  has  from  Eq. 
(5*1),  to  maximize 

A^5a  (6.1) 


subject  to 
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6 a TW  6 a < £ 2 . (6.2) 

a — 

Since  the  objective  i3  linear,  the  inequality  constraint 


can  be  replaced  by  an  equality.  Putting 


6a 


v^1 


2 

one  may  determine 

T -1 
A WaxA 


(6.3) 


(6.4) 


For  problems  that  are  quadratic  on  6a , one  has,  from 
Eq.  (5.1),  to  maximize 


AT6a  + ^6aT  AA6a  , (6.5) 


subject  to 


m p 

6a  W 6a  < C . 

Ol  — 


From  Eq.  (5.4),  one  finds 


-1, 


6a  * -(A I - 2vW  ) A. 


(6.6) 


(6.7) 


To  find  v for  this  problem,  an  iterative  scheme  is 
employed.  Let 


h(v ) = -(AA  - 2vW„ ) 1A. 


(6.8) 


Assume  a v°  to  get  the  direction  h=h(v°)  and  maximize  the 
objective  along  h to  find  6a°*Bh(v°),  where  B>0.  That  is, 
maximize 

6 g = BA^h  + |—  h^*AA*h,  (6.9) 

subject  to 

(6.10) 


B2hTW  h < e2. 
a — 


Now,  for  the  equality  case 

6 


B = 


/hVl  ■ 


(6.11) 
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For  the  Inequality  case, 

ATh  + 6 hT  AA  • h «=  0 


so 


(6.12) 


By  choosing  the  smaller,  nonnegative  6 from  the  above,  one 
has 


which  Is  the  Improved  value  of  v . Defining  the  Improved 
direction  h»h(\>),  the  procedure  Is  repeated  until  the 
change  of  v is  within  a prescribed  limit. 

For  the  POD  problem,  the  value  of  v chosen  in  this 
way  Is  only  an  approximation,  due  to  the  presence  of  the 
term , BT^fib,  in  Eq.  (4.7),  where  fib  is  still  unknown.  As 
noted  earlier,  in  the  presence  of  a -constraints , the  error 
from  the  optimal  would  be  increased.  This  latter  error, 
however,  would  not  be  significant  as  long  as  one  could 

still  move  on  the  tangent  plane  of  the  active  constraints. 

2 

Difficulties  can  be  encountered  when  £ is  too  large. 
Such  problems  in  the  one  dimensional  case  are  shown  in 
Figure  8.  If  £ is  small  enough,  these  difficulties  can  be 
avoided,  but  this  may  require  excessive  computing. 


0h(v°)  = -(AA  - 2vW  )“1A 

(X 


or  solving  for  v 


v * 


ft  hTAA»  h + hTA 

2fthTW  • h 
a 


(6.13) 


For  the  first  case  in  Figure  8,  there  seems  to  be  no 
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(1)  Missinp  completely 


(li)  Oscillating;  indefinitely 


S = starting  point 
N = next  point 


Figure  8 Difficulties  in  the  Choice  of  Stepsize 
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2 

solution  except  to  use  sufficiently  small  £ from  the 

2 

start.  For  the  second  case,  one  can  reduce  £ when 
oscillation  is  observed. 

It  is  evident  that,  when  oscillation  occurs  between 
two  points,  there  exists  a local  maximum  between  the  point 
points.  To  discover  this  oscillation,  the  inner  product 
• 6</ 1 V | I | | • | | 6a^  | | is  checked.  If  it 

is  larger  than  a prescribed  number,  say  0.9,  then  oscilla- 
tion is  clearly  occurring. 

As  an  example  of  these  difficulties,  one  may  have  a 
practical  situation  of  the  following  form  of  constraints 
(which  is  the  subproblem  of  POD) : 

a|z^(a)|  £Z2(a),  for  all  a c A.  (6.14) 

When  a^O,  there  is  no  difficulty.  When  a<0  however,  at  a 
point  where  z^(a)*0,  one  may  have  a local  maximum  point, 
where  the  derivatives  are  discontinuous  (not  satisfying 
the  assumptions  of  POD) . This  will  give  an  oscillatory 
behaviour  around  that  maximum  point . 

6.3*2  A Suggestion  on  the  Choice  of  y and  Wb 
In  the  gradient  projection  method  with  constraint 

error  compensation,  described  in  Section  2.4,  there  is 

2 

an  undetermined  constant  n in  the  stepsize  restriction 
given  by  Eq . (2.54).  From  Eq . (2.55),  it  is  easily  seen 
that  the  multiplier  y is  directly  related  to  the  stepsize. 

p 

Since  no  optimal  choice  of  n is  known,  the  precise 
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determination  of  y has  little  meaning  (Compare  with  the 
discussion  in  the  previous  section).  Hence  it  might  be 

p 

possible  to  make  y nonzero  by  assuming  that  n were 
adjusted  such  that  Sb^W^ 6b=n2 , still  satisfying  the 
condition  y(  <5bTWb 6b-n2)s0 , from  the  Kuhn-Tucker  theorem. 

In  reality,  one  can  think  that  the  stepsize  restriction  is 
always  active,  meaning  that  y should  be  positive.  Hence, 
in  the  following,  y>0  will  be  assumed.  Then,  from  Eqs. 


(2.5*0  and  (2.62), 


n2  > 6bTWb6b 


<5b^  Wh6b^  m 

5 + 5b2  Wh«b2,  (6.15) 

HY2 


since  4b^Wb6b2=0.  This  allows  an  explicit  lower  bound 
for  y such  that 


6blTWbfib2 

y > -7======-  » (6.16) 

2 / n2-6b2TWb6b2 

where  n2  > 6b2*^Wb6b2  Is  assumed.  Otherwise,  a compensa- 
tion procedure,  6b=6b2,  may  first  be  necessary. 

Prom  computational  experience,  Wb  is  usually  taken 
as  a diagonal  matrix  such  that 


^b^i  = (6.17) 

where  b^  denotes  the  magnitude  of  the  estimated  design 
component  and  Wb^  denotes  the  diagonal  element  of  Wb . Then, 
n is  taken  between  0.01  and  0.1,  at  the  beginning. 

Adjustment  of  y during  Iterations  often  leads  to  a better 


convergence . 
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6.3.3  Preliminary  Elimination  Procedure 

In  some  cases,  the  vectors  £*,  jele,  are  not  linearly 

J 

Independent,  giving  a singular  matrix  in  Eq.  (4.23). 

This  usually  happens  when  several  candidates  for  local 
maximum  points  are  assumed  at  the  beginning,  for  the  same 
g^ . With  poor  initial  estimates  for  designs  and  free  para- 
meters, there  can  be  too  many  constraints  violated.  To 
prevent  this  from  happening,  a preliminary  elimination 
procedure  may  be  applied. 

According  to  the  theorem  due  to  Fritz  John  (See  Sec- 
tion 2.4),  there  can  be  at  most  n("dimension  of  the  design 
variable  vector)  tight  constraints,  which  are  enough  to 
characterize  the  solution  point  b.  Therefore,  if  there 
were  more  constraints  violated  than  the  number  n,  it  might 

be  useful  to  consider  only  the  more  "significant"  con- 

T 6^"  J 

straints.  Suppose  g^,  is  violated  and  if  Vg^vf  ■ £*  £ >0, 

one  may  expect  that  after  an  iteration  this  violation 
would  be  compensated  for,  even  though  is  removed  from 
the  active  set.  Such  constraints  are,  thus,  deleted. 


6.3.4  Merge  Test 

A "merge  test"  is  given  as  often  as  required  to  check 
whether  the  assumed  local  candidate  maximum  points  for  g^ 
are  near  enough  to  give  a singular  matrix.  Since  one  has 
to  include  every  local  maximum  point  in  the  effort  to  find 
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the  global  maximum  point  for  g^ , there  is  a tendency  to 
have  more  inner  problems  than  are  really  necessary.  The 
merge  test  is  an  important  step  to  reduce  the  amount  of 
computation  by  eliminating  unnecessary  inner  problems,  as 
well  as  preventing  linear  dependency  of  the  gradients  of 
g^'s.  For  a check  of  nearness  of  two  vectors,  one  may  use 
the  maximum  norm  of  the  difference  vector. 

6.4  Integration  of  the  Algorithm 
6.4.1  Introduction 

Most  of  the  steps  described  above  need  to  be  repeated 
iteratively  to  obtain  accurate  results.  Consider  the  con- 
struction of  a sequence  b1  converging  to  the  limit  point  b 
by  using  another  sequence  6b^,  Sb^,***,  that  converges  to 
6b1.  One  then  puts  b1+3  = b1+6b1 . This  latter  construc- 
tion may  be  termed  a subalgorithm. 

Theoretically,  one  can  never  compute  the  exact  limit 
point.  In  practice,  therefore,  one  has  to  truncate  the 
construction  of  the  sequence.  If  one  lets  the  subsequence 
run  too  long  for  each  (outer)  iteration,  one  may  be  using 
more  computer  time  than  is  really  needed.  On  the  other 
hand,  if  one  truncates  too  soon,  he  may  degrade  conver- 
gence in  the  outer  iteration.  The  best  truncation  proce- 
dure is  not  known.  An  approach  to  introducing  a trunca- 
tion procedure  is  described  in  [36],  and  is  used  in  the 
present  work. 
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Since  one  can  measure  the  closeness  of  the  current 
iterate  to  the  limit  point  by  a scalar  quantity,  one  re- 
lates this  quantity  to  the  subalgorithm,  such  that  a cor- 
respondingly accurate  iterate  for  the  subsequence  is  ob- 
tained. In  reality,  this  quantity  is  not  known  until  the 
limit  point  is  obtained.  In  the  algorithm  model  shown  be- 
low, one  reduces  a scalar  quantity  from  a given  value  by 
a predetermined  ratio,  until  the  value  is  less  than  a pre- 
scribed small  number.  Such  an  algorithm  is  cited  in 
Figure  9 from  [36],  where  S(e,b*)  denotes  the  subalgorithm, 
c(b*)  is  a scalar  quantity  that  relates  the  closeness 
of  the  current  b*  to  b,  a>0,  0e(O,l),  eo>0  and  csc(0»6o) 
are  predetermined  numbers.  In  certain  cases,  the  regular 
reducing  of  e by  fixed  ratio  B would  result  in  slow  con- 
vergence. In  this  case,  the  algorithm  can  be  modified  by 
choosing  8 differently  after  each  iteration.  Such  a 
modification  is  Included  in  the  present  program. 

The  truncation  procedure,  as  explained  in  this  section 
using  a decreasing  e,  will  be  termed  an  e-procedure.  A 
convergence  property  for  this  algorithm  is  described  in 
the  same  reference,  under  assumptions  originally  due  to 
Zangwill  [48], 
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Figure  9 An  e -procedure  for  an  Algorithm  S(e , b*) 
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6. ^.2  Application  of  the  e -procedure  to  POD 

For  the  present  problem,  the  Newton's  method  employ- 
ed to  generate  b+fib  will  be  identified  as  a subalgorithm 
S(e,b),  with  the  e involved  in  the  stopping  criterion 
for  the  Newton's  method.  The  scalar  quantity  c(b*)  will 
be  identified  as  ||fib^|| . For  a normal  problem,  the  New- 
ton's method  should  generate  points  such  that  ||fib^||  is 
reduced.  But  for  POD,  several  other  procedures  are  involv- 
ed and  often  some  elements  b*  of  the  sequence,  generated 
during  the  iterations,  are  not  even  feasible  points,  ob- 
scuring the  convergence.  Convergence  is  further  obscured 
due  to  the  convergence  property  of  the  Newton's  method  it- 
self, especially  for  global  convergence. 

The  same  type  of  e -procedure  as  applied  to  the  New- 
ton’s method  above  can  be  used  in  the  truncation  of  the 
sequence  generated  by  the  maximization  procedure  and  in 
the  sign  check  of  the  various  multipliers  corresponding 
to  c -active  constraints.  For  the  latter  procedure,  the 
subalgorithm  is  similar  to  that  described  in  [36],  except 
that  for  the  present  computer  program  the  Lagrange  multi- 
pliers and  fib  are  calculated  by  another  subprocedure,  the 
Newton's  method.  Thus  the  overall  algorithm  for  POD  con- 
tains several  subprocedures,  interwoven  and  nested  in  such 
a way  that  it  is  much  more  complex  than  conventional  algo- 
rithms for  nonlinear  programming. 
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The  essential  part  of  the  flow  chart  of  the  POD  pro- 
gram Is  Shown  schematically  In  Figure  10.  Major  subalgo- 
rithms are  as  follows.  Loop  (1)  denotes  a Newton's  meth- 
od, whose  number  of  Iterations  depends  on  the  parameter 
c^,  which,  in  turn,  is  reduced  by  a factor  of  8c (0,1) 
when  b*  Is  near  the  solution.  Loop  (2)  denotes  the  sign 
check  of  various  mulitpliers.  Here,  also,  a slight  nega- 
tive value  of  multipliers,  by  an  amount  of  eA,  is  consid- 
ered to  satisfy  the  positivity  condition.  Whenever  the 
accuracy  of  the  maximum  point,  compared  with  that  of  design 
variable  is  poor,  the  maximization  procedure  denoted  by 
Loop  (3)  is  executed  until  an  accuracy  of  Cg  is  obtained. 

In  this  case,  the  branching  variable  MAX  is  set  equal  to 
one.  In  some  situations,  a regular  decrease  of  by  a 
constant  factor  B is  inefficient,  resulting  in  too  much 
time  in  Loop  (4),  which  decreases  c^  to  an  acceptable 
value  Eg.  To  check  whether  this  is  the  case,  a rapid 
decrease  in  eA  is  inserted  in  Loop  (*0  by  use  of  the 
branching  variable  IEADJ.  When  this  rapid  decrease  in  cA 
does  not  stop  the  program,  the  original  value  of  stored 
in  e ' is  recovered.  The  variables,  ITNEW,  INNER,  IMAX,  and 
ITER  register  the  number  of  iterations  for  the  subalgo- 
rithms, and  are  supposed  not  to  exceed  given  allowable 
numbers,  although  not  explicitly  expressed  in  the  flow 
chart . 


Figure  10 


Schematic  Flow  Chart  of  POD 


6.5  Conclusions 


A simplified  block  diagram  of  the  algorithm  developed 
in  Chapter  4 has  been  described.  Several  suggestions 
regarding  selection  of  parameters  in  the  algorithm  have 
been  given,  based  on  numerical  experience.  It  is  pointed 
out  that  the  conclusions  of  the  present  chapter  draw  on 
experimental  results.  Mathematical  Investigation  of  the 
algorithms  developed,  thus,  must  be  a subject  of  future 
study. 
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CHAPTER  7 

SUMMARY  AND  CONCLUSIONS 

In  this  thesis,  parametric  optimal  design  problems, 
motivated  from  the  concept  of  environmental  parameters,  have 
been  defined.  Solution  algorithms  have  been  developed  and 
applied  to  the  solution  of  a series  of  approximate  problems. 

The  basic  theoretical  tool  employed  is  the  implicit 
function  theorem,  applied  to  the  state  equations  and  the 
expansion  procedure  described  in  Sections  2.1.3  and  2.4.6. 
Algorithms  are  generated,  based  on  different  precision  of 
approximations.  In  Chapter  3,  a first  order  algorithm  is 
developed,  which  is  equivalent  to  the  alternating  maximi- 
zation and  minimization  technique  for  min-max  problems. 
Numerical  examples  show  that,  although  convergence  is  slow 
near  the  solution,  the  method  is  reliable. 

In  Chapter  4,  a second  order  approximation  has  been 
utilized.  The  inclusion  of  second  order  terms  is  moti- 
vated from  the  desire  to  adjust  the  maximum  points  in 
accordance  with  the  change  in  the  design  variable.  This 
latter  idea  leads  directly  to  sensitivity  analysis  for  the 
maximum  point.  It  has  been  pointed  out  that  this  analysis 
gives  a first  order  feedback  law  for  the  inner  problem. 
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This  relation  Is  essential  to  developing  various  other 
types  of  alternative  methods.  The  numerical  results  given 
In  Chapters  4 and  5 show  that,  near  the  solution,  this  pro- 
cedure converges  rapidly.  During,  the  sensitivity  analysis 
In  Section  4.2,  a compensating,  term  for  the  error  in  the 
nominal  maximum  point  is  found  to  be  the  same  as  that  ob- 
tained from  a second  order  method  of  solving  a nonlinear 
programming  problem.  This  method  of  maximization  gives  a 
different  point  of  view  for  the  damping,  factor,  as  appear- 
ed in  a damped  Newton’s  method  [45,47].  In  Chapter  5,  the 
possibility  of  some  other  variations  to  the  first  and 
second  order  algorithms  has  been  discussed,  with  numerical 
solutions  of  example  problems  presented. 

Numerical  experimentation  with  the  algorithms  shows 
that  their  success  is  problem  dependent.  That  is,  the 
structure  of  the  specific  problem  and  the  nature  of  the 
functions  involved  are  determining  factors.  Also,  it  has 
been  shown  In  Chapter  6 that  the  choice  of  parameters  in 
the  algorithm,  for  a g,iven  problem,  is  crucial. 

Table  17  summarizes  the  general  features  of  the  vari- 
ous algorithms,  based  on  the  numerical  experimentation. 

In  general,  the  algorithms  work  well,  although  there  is 
room  for  further  development.  It  is  hoped  that  the  present 
development  of  algorithms;  based  on  the  idea  of  expan- 
sion procedures,  sensitivity  analysis  with  error  compensa- 
tion, and  second  order  methods  of  maximization;  will 


motivate  theoretical  as  well  as  practical  investigation  of 
the  parametric  optimal  design  problem. 


Table  17  Comparison  of  Algorithms 


1-st  order  algorithm 

Hybrid  method 

2-nd  order  algorithm 

1.  reliability 

excellent 

excellent 

dependent  on  initial 
estimates  and  problem 
structure 

2.  convergence* 

slow 

good 

fast 

3.  programming 
difficulty 

simplest 

simple 

complicated 

computational 

effort 

small 

reasonable 

large 

5.  underlying 
methods 

gradient  projection 

2-nd  order  maximiza- 
tion, sensitivity 
analysis,  gradient 
projection 

2-nd  order  maximiza- 
tion, sensitivity 
analysis,  Newton's 
method 

6.  separate  maxi- 
mizations 

essential 

not  essential 

not  essential 

7.  overall  . 

justification* 

good 

better 

problem  dependent 

•The  rate  of  convergence  appears  to  approach  that  of  the  first  order  method 
(steepest  descent)  and  second  order  method  (Newton's  method)  of  NLP,  as  the 
nominal  design  approaches  to  the  solution.  This  is  because  the  solutions  of 
inner  problems  do  not  change  very  much  and  the  POD  problem  behaves  like  a NLP 
Droblem. 

+In  general,  quite  problem  dependent. 
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APPENDIX  D 


THE  VARIATION  OP  g IN  TERMS  OP  6 b AND  6a 


Since  it  is  assumed  that  the  system  of  state  equations 
3 In 

is  normal , i .e . , Yz  °>  2 “ z(b  ,a  ) , *g(z  ,b  ,a  ) « G(b  ,a  ) . 

It  is  necessary,  in  Section  4.2,  to  calculate  the  deriva- 
tives of  g(z(b  ,a  ) ,b  ,a  ) with  respect  to  b and  a.  By  use  of 
chain  rule  of  differentiation, 

= iS  + i£  ii 

3b  3b  3 z 3b 


3G 

3a 

32G 

3b7 

32G 
3 a 2 

32G 
3b  3a 

where 


i£  . iS  ii 

3a  + 3 z 3a 


iif  + 

3b7  + 

f4  + 


3ZT  3 2 


3^g  3 Z 3Z 

3b  3 Z 3b  + 3b  3z3b 


g 


3 2 Z 
3Z  3b2 


_Jl£  l2  + 

3a3z  3a 


3*g  3Z  3z 

+ 353z  ga  + ^ 


3ZT  3 2g 
3a  3z3a 

3zT  3 2 


3_Z  3^g  3Z  3g 

3b  3 Z2  3b  3Z 

3z^32g  3z  3j£  32Z 

3a  3 Z2  + 


3b  3a 


S 

3b  3z3a 


3a 
3 Z 


3Z  3a‘ 


3zT32g  3 Z 3g 

3b  3Z7  3a  + 3Z  3b 3a 


3ZZ 


3 2 Z 


(i£. 

v 3Z  3b 


32Z, 


U: 


_i£  v 

*zk  3bi3bj» 


etc . 


Then,  to  second  order  terms 
6 g * 6G 


3G 


3b 


b + 


3G 

3a 


32G, 


a + 


6 b^ 


3 2G , , 1 , 

JbdZ6a  + 26a  3a 


32G 


which  is  the  desired  expression.  Comparing  this  formula 
with  Eq.  (4.4),  the  matrices  B,  A,  SB”,  AA,  and  BA  are 
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identified . 

To  eliminate  various  derivatives  of  z in  the  most 


general  case,  differentiation  of  h(z,b,a)«0  yields 


ah  3z  ( ah  n 
az  ab  + ab  * u 


from  which  ^ is  obtained,  either  by  solving  sets  of  lin- 
ear equations  or  computing  if  the  inverse 

is  available.  Similarly, 

ah  az  . ah  _ az  /ah.-^ah 

az  3a  + 3a  °»  °r  da  '3Z}  3a* 


Differentiating  again. 


, azTa2h  a2hv az  ah  a2z 

' ab  az2  abaz'ab  az  ab2 

Hence , 

a2h  / a2h  az  azT  a2hv  a zT a 2 h az 
ab7"  “ ^abaz  ab  ab  azab;  “ ab  Tz7  ab 


aji  a2z 
az  ab2 


a2h  azT  a2h  n 
ab7*  + ab  azab  “ u 


Similarly , 


ah  a2z 
az  a a7 


/ a2h  az  azTa2h  v 
V dodz  3a  3a  3z3a; 


a zTa 2h  az 
a q a z 2 aa’ 


and 


ah  a2z 
az  abao 


a2h  / a2h  az  +azTa2h  x 
abaa  “'•3b3z  3a  3b  3z3a7 


azT a2h  az 
ab  az*  3a 


It  is  observed  that  if  the  dimension  of  the  state  variable 


z is  small,  it  may  be  computationally  more  profitable  to 
9 h 1 

compute  (y^)~  . In  general,  however,  the  computation  is 
very  time-consuming.  As  a close  approximation,  if  the 
state  equations  are  linearized  as  follows, 
ah,  , ah, . , ah,  a 

+ j^6a  - 0, 


then,  since 
computation 


the  second  derivatives  of  z are  zero,  the 
is  enormously  reduced. 
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APPENDIX  E 

EQUIVALENCE  OF  LAGRANGE  MULTIPLIERS  OF  NLP  AND 
THOSE  OF  ITS  FIRST  ORDER  APPROXIMATE  PROBLEM 

The  Lagranpe  multipliers  for  the  approximate  problem 
are  given  by  Eq.(2.56).  At  the  solution  b,  Ag  ■ 0.  Hence, 


The  subscript  ♦ denotes  the  active  constraints  at  b. 

Now  consider  the  original  NLP  given  by  Eqs.  (2.43)  and 
(2.44)  without  equality  constraints.  At  the  solution  b, 
Kuhn-Tucker  theorem  assures  the  existence  of  multipliers 
such  that 

If  ’ ° * H ♦ >TH.  (E-2> 

where  H = f + xTg>  and  the  only  possibly  nonzero  elements 
are  kept  for  X.  That  is,  if  g^<0,  X^=0  so  the  nonzero  ele- 
ments of  X correspond  to  the  active  constraints  only.  In 
the  notation  of  and  calculated  at  b, 


Since  the  existence  of  X is  insured,  postmultiplying 


x = -m^1  i*Tvr1iJ, 

where  and  and  l* 


are  calculated  at  b. 


(E.l) 


(E.3) 


through  by  W”^i^  and  transposing,  one  has 


(E . 4 ) 


which  gives  multipliers  X,  same  as  obtained  in  Eq.  (E.l). 
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